Let's imagine that we have a straight-line fitting function
         y = mx + b   where m is the slope and b is the intercept.

m = -1./2.; b=8.;σ=0.3;

Let's imagine that we have 31 data points where the x coordinate ranges from -10 to 20.

num=31;

x = Table[i,{i,-10,num}]

Let's randomly smear our data in the "y" direction with an 
uncertainty of σ = 0.30.

rarray=RandomVariate[NormalDistribution[0,σ],num];
Histogram[rarray,10]

Okay . . . .   let's smear the data in the "y" direction.

y = Table[m x[[i]]+b + rarray[[i]],{i,1,num}];
data=Table[{x[[i]],y[[i]]},{i,1,num}];

Okay . . . . let's do a least-squares fit and find the slope "d" and the intercept "c."

f1 = FindFit[data,c+d t,{c,d},t]
Plot[c + d t/.f1,{t,-12,num},Epilog->Prepend[Point/@data,PointSize[0.013]],AxesOrigin->{0.,0.}]

dev = Table[y[[i]]-(c+d x[[i]] /.f1),{i,1,num}];

Let's plot the deviations of the data point from the least-squares fit to a straight line.

Histogram[dev,10]

Let's calculate the (χ^2).

chisqmin=(∑_(i=1)^num)(((((dev[[i]])/σ)))^2)

Let's use the traditional least-squares technique to find the:
intercept:  a ± δa    and
slope:         b ± δb.

delta =(∑_(i=1)^num)(1/(σ^2)) (∑_(i=1)^num)(((x[[i]])^2)/(σ^2)) - ((((∑_(i=1)^num)((x[[i]])/(σ^2))))^2);
a1 = (1/delta)((∑_(i=1)^num)(((x[[i]])^2)/(σ^2)) (∑_(i=1)^num)((y[[i]])/(σ^2)) - (∑_(i=1)^num)((x[[i]])/(σ^2)) (∑_(i=1)^num)((x[[i]] y[[i]])/(σ^2)));
b1 = (1/delta)((∑_(i=1)^num)(1/(σ^2)) (∑_(i=1)^num)((x[[i]] y[[i]])/(σ^2)) - (∑_(i=1)^num)((x[[i]])/(σ^2)) (∑_(i=1)^num)((y[[i]])/(σ^2)));
deltaa = (Sqrt[(1/delta)(∑_(i=1)^num)(((x[[i]])^2)/(σ^2))]);
deltab = (Sqrt[(1/delta)(∑_(i=1)^num)(1/(σ^2))]);
Print["a =  ",a1, "  ± ",deltaa]
Print["b = ",b1, " ± ",deltab]

Lo . . . and behold . . . .  it's the same results for the slope and intercept that we obtained using FindFit.  No surprises here.

Now, let's generate the (χ^2)'s as a function of a and b, and plot them.

chisqA[ax_]:= (∑_(i=1)^num)((((y[[i]] - (ax+b1*x[[i]])))^2)/(σ^2))
chisqB[bx_]:= (∑_(i=1)^num)((((y[[i]] - (a1+bx*x[[i]])))^2)/(σ^2))
Plot[chisqA[ax],{ax,a1-deltaa,a1+deltaa},AxesOrigin->{a1,chisqmin}]
Plot[chisqB[bx],{bx,b1-deltab,b1+deltab},AxesOrigin->{b1,chisqmin}]

Notice that as we move from (Subscript[a, min ])→(Subscript[a, min ])±δa, the (χ^2) changes from (Subscript[χ, min]^2)→(Subscript[χ, min]^2)+1.  Likewise, when we move from (Subscript[b, min ])→(Subscript[b, min ])±δb, the (χ^2) changes from (Subscript[χ, min]^2)→(Subscript[χ, min]^2)+1.

Given the above curves, the uncertainties δa and δb can be determined by solving for the intersection of (Subscript[χ, min]^2)+1=(χ^2)((Subscript[a, min])+δa) and the same for (Subscript[χ, min]^2)+1=(χ^2)((Subscript[b, min])+δb)

amin=c/.f1;
bmin=d/.f1;
fa = FindRoot[chisqA[ax]==chisqmin+1.,{ax,amin}];
errora=Abs[(ax/.fa )- amin];
fb = FindRoot[chisqB[bx]==chisqmin+1.,{bx,bmin}];
errorb=Abs[(bx/.fb )- bmin];
Print["a = ",amin  ," +/- ", errora];Print["b = ",bmin  ," +/- ", errorb];

Let's generate (χ^2)'s from multiple experiments and let's see how they are distributed.

chisqArray={};
m = -1./2.; b=8.;σ=0.3;
num=31;
x = Table[i,{i,-10,num}];
Do[{rarray=RandomVariate[NormalDistribution[0,σ],num];
y = Table[m x[[i]]+b + rarray[[i]],{i,1,num}];
data=Table[{x[[i]],y[[i]]},{i,1,num}];
f1 = FindFit[data,c+d t,{c,d},t];
dev = Table[y[[i]]-(c+d x[[i]] /.f1),{i,1,num}];
chisq=(∑_(i=1)^num)(((((dev[[i]])/σ)))^2);
 chisqArray=Append[chisqArray,chisq]},{ii,1,100}]
g1=Histogram[{chisqArray,10},20,PlotRange->{{0,60},All},AxesOrigin->{0,0}]

probChi[x_,ν_]:=(((x^(((1/2)(ν-2)))) (E^(-(x/2))))/((2^(ν/2))  Gamma[(ν/2)]));

g2=Plot[220*PDF[ChiSquareDistribution[29],x],{x,0,60},PlotRange->All]

g3=DiscretePlot[220 PDF[ChiSquareDistribution[29],x],{x,1,60,2.},PlotRange->All]

Show[g1,g3]

At what value of (χ^2) does the area underneath the Probability((χ^2)) approach 95%?

(∫_42.556^∞)probChi[chi2,29]\[DifferentialD]chi2//N

This occurs at   (χ^2)=42.556   for 29 degrees of freedom.